import open3d as o3d
import laspy
import numpy as np
import matplotlib.pyplot as plt

# ----------------------------
# 1. 读取 LAS 文件
# ----------------------------
las_path = "/mnt/d/temp_files/part_pointCloud.las"  # 请根据实际情况修改路径
las = laspy.read(las_path)

points = np.vstack((las.x, las.y, las.z)).T  # shape: (N, 3)

print(f"原始点数: {points.shape[0]}")

# ----------------------------
# 2. 创建 Open3D 点云并降采样
# ----------------------------
pcd = o3d.geometry.PointCloud()
pcd.points = o3d.utility.Vector3dVector(points)

# 设置体素大小（单位：与 LAS 坐标一致，例如米）
voxel_size = 0.5  # 可根据点云密度调整，如 0.2, 1.0 等

pcd_down = pcd.voxel_down_sample(voxel_size=voxel_size)
downsampled_points = np.asarray(pcd_down.points)

print(f"降采样后点数: {downsampled_points.shape[0]}")

# ----------------------------
# 3. 绘制 X, Y, Z 方向的直方图
# ----------------------------
fig, axes = plt.subplots(1, 3, figsize=(18, 5))

# X 方向
axes[0].hist(downsampled_points[:, 0], bins=100, color='red', alpha=0.7)
axes[0].set_title('X Coordinate Distribution')
axes[0].set_xlabel('X')
axes[0].set_ylabel('Frequency')

# Y 方向
axes[1].hist(downsampled_points[:, 1], bins=100, color='green', alpha=0.7)
axes[1].set_title('Y Coordinate Distribution')
axes[1].set_xlabel('Y')
axes[1].set_ylabel('Frequency')

# Z 方向（高度）
axes[2].hist(downsampled_points[:, 2], bins=100, color='blue', alpha=0.7)
axes[2].set_title('Z Coordinate (Height) Distribution')
axes[2].set_xlabel('Z')
axes[2].set_ylabel('Frequency')

plt.tight_layout()

# 可选：保存图像
# plt.savefig("xyz_histograms.png", dpi=150)

# 显示图像
plt.show()